Self-consistent theory of pair distribution functions and effective interactions in 

quantum Coulomb liquids 
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We use a density-functional theoretical approach to set up a computationally simple self-consistent 
scheme to calculate the pair distribution functions and the effective interactions in quantum Coulomb 
liquids. We demonstrate the accuracy of the approach for different statistics and space dimension- 
alities by reporting results for a two-dimensional electron gas and for a three-dimensional boson 
plasma over physically relevant ranges of coupling strength, in comparison with Monte Carlo data. 
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Exchange and correlations in systems of interacting 
electrons have been a focus of interest in many-body 
physics for many decades. An important manifestation 
of these effects is the equilibrium distribution of elec- 
tron pairs, which provides a picture of the short-range 
order in the system. Knowledge of this distribution for 
a homogeneous electron fluid as a function of its density 
determines its ground-state energy and is essential in the 
construction of energy functionals for the study of atoms, 
molecules and solids in applications of density functional 
theory (DFT) transcending a local-density approxima- 
tion p]. Recent developments in DFT have drawn at- 
tention to the utilization of inhomogeneous electron-pair 
densities for such studies || . There also is an interest for 
fluids of charged bosons in quantum statistical mechan- 
ics, e.g. in regard to condensates ofpoint-like Cooper 
pairs as models for superfluid states g . 

We present in this Letter a novel theory of the pair dis- 
tribution function g(r) in homogeneous fluids of charged, 
point-like fermions or bosons. We start from some ba- 
sic results of DFT [Q and tackle the evaluation of g{r) 
through the solution of two-particle scattering problems 
which are governed by effective scattering potentials em- 
bodying the many-body effects. This viewpoint, which 
was first taken by Overhauser || and further developed 
in subsequent studies @, [J, is here implemented through 
a self-consistent inclusion of exchange and correlations 
leading back, when a non-self-consistent linear-response 
approximation is made, to the spin-dependent effective 
interactions derived in the early work of Kukkonen and 
Overhauser ||. Calculations on fluids of electrons in di- 
mensionality D = 2 and of charged bosons in D — 3 are 
shown to yield results in excellent agreement with avail- 
able Quantum Monte Carlo data over a physically sig- 
nificant range of coupling strength. As previously noted, 
the method can be extended to pair distributions in in- 
homogeneous Coulomb fluids |Q. 

We consider a quantum fluid of point-like particles hav- 
ing charge e and mass m at zero temperature, neutral- 
ized by a uniform charged background. The fluid con- 
sists of different components at average density n a (two 
spin components for spin-1/2 fermions, for instance). As 
already stated, our aim is to use DFT for building a self- 



consistent theory that gives the pair distribution func- 
tions as output. From the Hohenberg-Kohn theorem [0, 
the ground-state energy of the fluid in the presence of a 
set of external potentials V„ xt (r) can be written as 

E gs [Mr)}} =T s [{n ff (r)}] + E H [{n a (r)}} 

+ W ' d D v C'(>-)An a (r) + E Q [{n a (r)}} (1) 

a J 

where {n a (r)} is the set of densities of all components, 
An Q (r) = n a (r) — n a , and T s is the ideal kinetic energy 
functional. The Hartree term E n is given by 
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J d D r J d D r' v(\r - r'\) An a (r)An f3 (r') 
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where v(\r — r'|) = e 2 /|r — r'|. The last term in Eq. 
is the exchange-correlation energy functional, containing 
the quantum many-body (QMB) effects. In Eqs. (|l|) and 
(0) the presence of a neutralizing background has been 
taken into account. 

The quantity n a [g ai {r) — 1], where g a -y(r) is the 
component-resolved pair distribution function, can be 
viewed as the distortion that a particle of the fluid (of 
7-type at position r = 0) produces in the density pro- 
files Here, g a j(r) is defined through the average 
number of cn-type particles inside a spherical shell of ra- 
dius r and thickness dr centered on a 7-type particle lo- 
cated at the origin, which is given by n a g a -y( r )^D rD1 dr 
with CI2 = 27r and ^3 = An. The appropriate ground- 
state energy functional is obtained from Eq. (|lj) with 
the formal replacements V^ xt (r) — * v(r) and An a (r) — > 
n a [g a j( r ) — 1]- Finally, the QMB energy functional can 
be written via an adiabatic connection formula, 

2 
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x n a (r)n p {r') [gK (r,r') - l] , 



(3) 



where n a (r) — n a g a ^ (r) and fl„ J g 7 (r,r / ) measures the 
probability of finding two particles with indices a and f3 
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at r and r' when a 7-type particle is at the origin, the 
interaction potential being vx (| r — r' |) = Ae 2 /|r — r'|. Of 
course, g^^r, r') depends functionally on {n a (r)}. 

The Kohn-Sham mapping Q ensures that n a (r) can be 
built from Kohn-Sham scattering orbitals ^^(r) which 
satisfy the following set of Schrodinger equations: 



-fv^(r) 



(4) 



Here, r is the relative distance of two particles, /1 = m/2 
is the reduced mass, and £k = fi 2 k 2 /2/x with tik the rel- 
ative momentum. The scattering potential V^(r) in 
Eq. (|J) can be obtained from the first functional deriva- 
tive of E gs — T s with respect to n a (r) , 



VS{r) = v(r) d D r'v(\r - r'|)A^(r') 



SE, 



(7) 



Sn a (r) ' 
(5) 



Finally, the pair distribution functions can be obtained 
from the Kohn-Sham scattering states as 
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where the sum runs over all occupied states labelled by 
k. The occupation factors T^ 7 depend on the statistics of 
the fluid (see the discussion below). Notice that Eq. (|^) 
guarantees positive dcfinitcncss of g a -y{r). 

In the above formal development, the functional de- 
pendence of the QMB energy on density is not known and 
we have to resort at this point to some approximations. 
Their goodness can only be gauged a posteriori. Firstly, 
the function (^^(r, r') in Eq. (||) involves three-body 
correlations and would lead us into a hierarchy of higher- 
order correlation functions. The simplest way of truncat- 
ing this hierarchy is to replace ff^O, r') by fl^Gr-r'l), 
in analogy with what has been done in treating the equa- 
tion of motion for the Wigner distribution function in the 
presence of external potentials [ flO| . Secondly, we expand 
the QMB energy in a functional Taylor series in powers 
of An a (r) up to second order terms. With the definition 



F/3(|r-r'|)= S 2 E Q [{n a (r)}] 
6n a (r)8np(r') 



(J) 



(An„(r)(=0 



we find in Fourier transform 



V^(q) = v(q) 



G af3 (q)}An (q) . (8) 



Here, G a p(q) = —f al3 (q)/v(q) are the so-called local field 
factors, defined in term of the Fourier transform f a @(q) 
of the QMB kernels in Eq. (0), and v(q) is the Fourier 
transform of the Coulomb potential (i.e. v(q) — 4ire 2 /q 2 
in D = 3, v(q) = 2-Ke 2 /q in D = 2). It may be worth 
exploring in the future alternative approximations to a 
truncated expansion of the QMB energy. 



The approximations that have led us to Eq. (||) can be 
justified in a weak-coupling regime and indeed Eq. (||) 
yields back the effective electron-electron interactions of 
Kukkonen and Overhauser pi when Anp(q) is related to 
the scattering potential V£s (<?) by linear response theory. 
We propose instead to carry out a self-consistent calcula- 
tion of the pair distribution functions and of the effective 
interactions through the solution of the set of equations 
©' © and @- For this purpose the quantities Anp(q) 
in Eq. (ph should be written in the form 



Ana(q) = (npln^ 12 [Sp 7 (q) - 8 01 ] 



(9) 



where Sp^^q) are the partial structure factors, related to 
the pair functions by 

s i3t(q) = <Vy + V n /3 n 7 J d ° r tehM - 1] ex P (-«q ■ r ) • 

(10) 

The uselfulness of such a self-consistent approach will 
be explored in the numerical calculations reported fur- 
ther below. Here we remark that setting G a p{q) = in 
Eq. (0) gives back the self-consistent Hartree approxima- 
tion (HA) , that we have found to yield quite satisfactory 
results for g(r) in the 3D electron gas for values of the 
coupling strength up to at least r s = 10 0. 

Let us examine the application of the above self- 
consistent scheme (SCS) to a 2D paramagnetic electron 
gas (EG), where correlations are stronger than in 3D. 
The Greek indices become spin indices taking the val- 
ues a = ±1. For the details of the summation pro- 
cedure in Eq. (^J), the reader is referred to Refs. || 
and 0]. The necessary input are the local- field fac- 
tors G alT '(q), for which we use Quantum Monte Carlo 
(QMC) dat a [pd[ as described by interpolation formu- 
lae in Ref. ]12[. These refer to the charge-charge and 
spin-spin field factors, G + (q) — [G^(q) + G]i(q)]/2 
and G-(q) = [G-\-\(q) — G-\i(q)]/2, which in the long- 
wavelength limit satisfy the compressibility and suscep- 
tibility sum rules, G+(q) — > (kq 1 — KT r )/[n 2 v(q)] and 

G-(q) -> ^Kxp 1 - Xs r )/ V (q)- Here > K o and \p are 
the compressibility and the Pauli susceptibility of the 
ideal Fermi gas, k and \s are the corresponding quan- 
tities for the EG, and /j,b is the Bohr magneton. We 
recall that these expressions, though strictly valid only 
in the thermodynamic limit, give in practice a good ac- 
count of G+(q) over a range of q extending almost up to 
2k F @ g|. 

FigureTjTreports our SCS results for the spin-averaged 
pair distribution function g(r) = [g]-](r) + g^i(r)]/2 in 
the 2D EG at r s = 5 and 10, where r s = (nna^)^ 1 ^ 2 is 
the usual coupling-strength parameter with as the Bohr 
radius. The results of a QMC study by S. Moroni (pri- 
vate communication) and of our previous HA approach 
are also shown in Figure |l|. It is seen that inclusion of ex- 
change and correlation in the present SCS reproduces the 
formation of a first-neighbor shell with increasing cou- 
pling strength, which is missed in the HA |t|, and yields 
satisfactory quantitative agreement with the QMC data. 
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Figure |2] reports the SCS results for the spin-symmetric 
component of the effective electron-electron interaction, 
Vks(r) = [Vjgir) + V^{r)}/2, for the 2D EG at r s = 5, 
as well as its Fourier transform Vks(<?) (shown in the 
inset). We find very significant changes from the HA 
and indeed the original treatment of exchange and corre- 
lations in a linear-response approximation by Kukkonen 
and Overhauser [^| are not very far from our SCS results. 

We should remark at this point that the present SCS 
is not as accurate in reproducing the spin-resolved pair 
distribution functions of the 2D EG and hence the spin- 
spin effective interaction. Numerical evidence that an 
empirical inclusion of higher-order terms in the scatter- 
ing potentials can yield full agreement with the QMC 
data, as well as extensions to larger values of r s , will be 
reported in future work. 

Here we discuss instead how our approach can be ex- 
tended into a fully self-consistent scheme (FSCS), in 
which the local field factors are self-consistently deter- 
mined over the relevant g-range during the calculation 
rather than taken as input from QMC. We need for this 
purpose a closure relation between G a p{q) and S a /3(q), 
and the crucial point is that it should self-consistently 
satisfy the thermodynamic (compressibility and suscep- 
tibility) sum rules. Such a closure can be obtained by 
using Eq. (||) in Eq. (^), with the aforementioned ap- 
proximation g^ ( g 7 (r, r') ~ 5q^(|i" — r'|). The details of 
this calculation will be reported elsewhere. The final re- 
sult is 

G a p(q) = D af3 g a p{q) . (11) 
Here, the differential operator D a p is defined by 

D a p = l + n a - \-np- V-n a np- — , (12) 

on a onp 2 on a onp 

while G a (3{q) is given by 

r <„\- 1 1 f 2 ^ f V W) 
G a p(q) = — / dX — — n — — 

e z J Q J (2tt) u v{q) 
x [5^(|q + q'|)-M , (13) 

with S*p(q) being the partial structure factor at coupling 
constant A. Although these expressions are strictly cor- 
rect only in the long- wavelength limit, they yield a good 
account of QMC data on the local field factors over the 
relevant range of q jn], [l2| , as already noted. 

It can be seen from Eqs. (|TT|)-([T3|) that, in order to 
satisfy the susceptibility sum rule in the paramagnetic 
EG, one needs to move into a partially spin-polarized 
state and indeed evaluate the full range of spin polar- 
ization up to the ferromagnetic state. For a simple il- 
lustration of our method we have therefore considered 
a 3D fluid of spinless charged bosons (CBF), where the 
proposed FSCS simplifies drastically. At zero temper- 
ature all bosons in the reference Kohn-Sham ideal gas 



are in the k — state, so that g(r) cx |3>k=o(r)| 2 . 
The scattering orbital $k=o(r) is a spherically symmet- 
ric function and the Kohn-Sham Schrodingcr equation 
becomes equivalent to the Euler-Lagrange equation for 
g(r) as obtained from the variational principle using 
the von Weizsacker functional for T s |13[, T s [g(r)] = 
(h 2 n/8/j,) J d D r\Vg(r)\ 2 /g(r). Furthermore, the local 
field factor G[q) from Eqs. (pT|)-(^3|) satisfies the com- 
pressibility sum rule, i.e. \\m q ^QG{q) = — [n 2 Kv(q)]~ 1 
with k being the compressibility of the interacting Bosc 
liquid. 

We have solved the FSCS based on Eqs. ({§), (§), 
(§) and ©-© for a 3D CBF at coupling strength 
r s = (47rna^/3) -1 / 3 up to 20. The main results are 
shown in Figures | and f§ In Fi gure [| we compare the 
FSCS g(r) with QMC data by Moroni et al. @ and with 
SCS results where the QMC data for the local field fac- 
tor [ fL4| have been used as input. Our results for g(r) are 
in excellent agreement with QMC and with each other, 
the implication being that our self-consistent determina- 
tion of G(q) from the compressibility sum rule also ac- 
counts with sufficient accuracy for this function over the 
relevant range of q (see also the inset in Figure |i| where 
the FSCS G(q) is compared with the QMC data). In 
the inset in Figure § we include, with the FSCS results 
and the QMC data for g(r) at r s = 20, also the results 
obtained in the HA and those of Apaja et al. || based 
on a hypernetted chain (HNC) approximation. Again, 
the HA is not able to reproduce quantitatively the for- 
mation of a first-neighbor peak with increasing coupling. 
Finally, in the main body of Figure || we illustrate the 
self-consistency that we have obtained in the values of 
the compressibility k from the ground-state energy and 
from G(q), as well as their excellent agreement with QMC 
data by Moroni et al. Q . 

In summary, we have proposed a self-consistent ap- 
proach by which the pair correlations (and related quan- 
tities such as the internal energy) of a quantum Coulomb 
fluid can be determined, ultimately using the partial 
mean densities of its components as the only input. 
We have examined its usefulness in different statistics 
(fermions and bosons) and space dimensionalities (D = 2 
and 3) in comparison with available QMC data, extend- 
ing up to relatively large values of the Coulomb coupling 
strength. We have also highlighted the role of exchange 
and correlations in determining the emergence of liquid- 
like structure with increasing coupling strength through 
the formation of a first-neighbor shell and further oscil- 
lations in the pair distribution function. 

This work was partially supported by MIUR through 
the PRIN2001 program. We are grateful to Prof. G. 
Vignale for stimulating discussions and to Dr. S. Moroni 
for giving us access to his QMC results for the 2D EG 
prior to publication. 
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FIG. 1: The pair distribution function g(r) in a 2D EG at r s = 5 and 10, as a function of r/(r s as). The results of the SCS (full 
line) and of the HA (dotted line) are compared with QMC data (crosses). The curves at r s = 10 have been shifted upwards by 
0.4. 



6 




_0 3 I i i i i i i I 

" 1 2 3 4 5 6 7 

r /(r s a B ) 



FIG. 2: The spin-symmetric effective electron-electron interaction Vks(r) in a 2D EG at r a = 5, as a function of r/(r s as). The 
result of the SCS (full line) is compared with the HA (dotted line) and with the Kukkonen-Overhauser result (dashed line). 
The inset shows the Fourier transform of Vks(r) as a function of qr a aB- The black dot gives the long- wavelength limit of the 
Thomas-Fermi theory for the electron-test charge interaction. 
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FIG. 3: The pair distribution function g(r) in a 3D CBF at r s = 10 and 20, as a function of r/(r s os). The results of the 
FSCS (full line) and of the SCS (dashed line) are compared with QMC data (crosses). The curves at r s — 20 have been shifted 
upwards by 0.4. In the inset, the HA result (dash-dotted line) is compared with the FSCS (full line), the QMC data (crosses), 
and the HNC of Apaja et al. (dotted line). 
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FIG. 4: (nft) -1 in Ryd as a function of r s for a 3D CBF. The results of the FSCS from the ground-state energy (full line) and 
from the self-consistent G(q) (empty boxes) are compared with QMC data (triangles). The inset shows G(q) as a function of 
qr a dB at r s — 10: the FSCS (full line) is compared with QMC data (dashed line). 



